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Abstract 

This  work  presents  a  macro-level  model  to  determine  the  performance  characteristics  of  solid  oxide  fuel  cells  (SOFCs).  Activation,  ohmic, 
and  concentration  polarizations  are  considered  the  main  sources  of  irreversibility.  The  Butler- Volmer  equation,  the  dusty  gas  model,  and  Ohm’s 
law  were  used  to  determine  the  polarization  terms.  Tafel  equation,  linear  current  potential,  Fick’s  model,  and  Stefan-Maxwell  model  were 
quantitatively  analyzed  as  well.  Performance  curves  were  calculated  for  hydrogen,  methane,  and  carbon  monoxide  as  pure  fuels  at  different 
conditions.  A  surface  plot  for  each  polarization  term  allowed  us  to  analyze  the  contribution  to  voltage  loss  as  a  function  of  temperature 
and  current  density.  The  calculations  presented  in  this  paper  involved  the  creation  of  many  computational  tools.  One  sample  code  of  these 
algorithms  is  included,  but  all  algorithms  used  in  the  paper  are  available  from  the  authors. 

©  2004  Elsevier  B  .V.  All  rights  reserved. 
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1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  operate  at  high  temperature 
(800-1200  °C),  which  makes  possible  the  use  of  a  variety  of 
fuels  (hydrocarbons),  cogeneration,  and  bottoming  cycles. 
The  Energy  &  Environmental  Research  Center  (EERC)  at 
the  University  of  North  Dakota  (UND)  designed  a  thermally 
integrated  biomass  SOFC-gasification  system  (BG-SOFC). 
The  producer  gas  (a  mixture  of  hydrogen,  methane,  carbon 
monoxide,  carbon  dioxide,  small  amounts  of  tars  <1%,  and 
steam)  from  the  gasifier  will  be  fed  into  the  SOFC,  and  the  hot 
effluent  gases  from  the  anode  and  cathode  will  be  reinjected 
into  the  gasifier  (at  different  gasification  zones),  maintain¬ 
ing  the  reactor  temperature.  At  high  temperatures  (^600  °C), 
methane  will  be  internally  reformed  in  the  SOFC  and  car¬ 
bon  monoxide  will  be  converted  to  hydrogen  (assuming  that 


*  Corresponding  author.  Tel.:  +1  7017779495. 

E-mail  address:  ehernand@und.nodak.edu  (E.  Hernandez-Pacheco). 

0378-7753/$  -  see  front  matter  ©  2004  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 1016/j.jpowsour.2004.06.05 1 


the  water-gas  shift  reaction  will  always  be  at  equilibrium) 
[1,2].  Conversion  efficiencies  are  expected  up  to  45%  based 
on  modeling  analyses  [3].  The  accurate  prediction  of  SOFC 
operating  conditions  is  important  for  attaining  high  efficien¬ 
cies  and  successfully  integrating  both  systems  (BG-SOFC). 
The  spectrum  of  possible  operating  conditions  can  be  deter¬ 
mined  based  on  the  performance  characteristics  of  the  fuel 
cell.  These  characteristics  will  allow  us  to  determine  a  set  of 
operating  conditions  (temperature,  pressure,  inlet  composi¬ 
tions)  for  optimal  performance. 

The  fuel  cell  performance  is  subject  to  the  second  law  of 
thermodynamics  so  that  losses  in  the  system  are  inevitable. 
For  a  fuel  cell,  the  main  source  of  irreversibility  comes  from 
(1)  ionic  resistance  through  the  electrolyte  and  electronic  re¬ 
sistance  on  the  electrodes  and  interconnects,  (2)  activation 
energy,  and  (3)  the  diffusion  of  gases  to  the  reaction  sites. 
These  losses  can  be  determined  numerically  or  experimen¬ 
tally  (see  [4,5]  for  experimental-related  literature).  However, 
numerical  models  are  preferred  for  economical  reasons  and 
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because  sensitivity  to  different  parameters  can  be  investigated 
with  the  same  model.  There  are,  in  general,  three  types  of  nu¬ 
merical  models  used  for  determining  the  polarization  terms  in 
SOFCs:  micro-level  models  [6-10],  macro-level  models  [11- 
14],  and  semiempirical  models  [15-17].  Results  presented  in 
this  paper  are  based  on  a  macro-level  model.  This  model 
is  based  on  the  assumption  that  the  reactions  occur  at  the 
electrode-electrolyte  interface,  and  therefore  no  analysis  at 
the  microscopic  level  is  required.  This  assumption  is  valid 
for  pure  electronic  conductors  but  is  unrealistic  for  cermet- 
type  anodes,  because  of  the  intermixing  in  the  electrode¬ 
electrolyte  interface.  However,  based  on  the  cermet-modeling 
results  by  Chan  and  Xia  [1 8] ,  it  can  be  argued  that  most  of  the 
polarization  contribution  occurs  at  the  interface,  especially 
for  anode-supported  cells  with  thick  anodes  (^600  pan).  The 
objective  of  this  paper  is  to  present  a  set  of  computational 
tools  for  predicting  the  polarization  terms  in  anode- supported 
cells  using  a  macro-level  model.  The  model  presented  here 
can  be  extended  into  a  system  model  that  includes  tempera¬ 
ture  and  flow  distributions,  and  this  work  will  be  presented 
in  the  future. 

Section  2  describes  different  approaches  to  calculate 
the  polarization  terms  (i.e.,  activation,  concentration,  and 
ohmic).  Numerical  results  from  the  Tafel  equation  and  lin¬ 
ear  current  potential  were  compared  with  the  Butler- Volmer 
equation.  Semiempirical  correlations  based  on  experimen¬ 
tal  kinetic  data  were  also  compared  to  the  Butler- Volmer 
equation.  Fick’s  model,  the  dusty  gas  model,  and  the  Stefan- 
Maxwell  model  were  used  to  determine  the  diffusion  of  re¬ 
actants  into  the  porous  material  (concentration  polarization). 
For  the  ohmic  polarization,  different  expressions  of  the  elec¬ 
trolyte’s  conductivity  were  analyzed.  Section  3  discusses  the 
results  obtained  from  the  macro-level  model.  The  perfor¬ 
mance  characteristics  were  calculated  at  different  operat¬ 
ing  conditions  and  for  three  electrochemically  active  fuels 
(i.e.,  hydrogen,  carbon  monoxide  and  methane).  Appendix  A 
describes  the  formulae  required  to  calculate  concentration 
polarization.  A  set  of  rational  approximations  for  the  most 
commonly  used  binary  systems  in  SOFC  technology  was  in¬ 
cluded.  A  sample  code  implementation  of  the  dusty  gas  model 
using  Mathematica®  is  presented  in  Appendix  B. 


2.  A  review  of  polarization  in  SOFCs 


given  by 

dG  =  -SdT+VdP 

+  (- VA/XA  -  vB/xB  +  vc/xc  +  vdMd)  (2) 


dG  =  -SdT+  VdP  +  dGZp  (3) 

where  /x  represents  the  chemical  potential,  and  d£  is  a  propor¬ 
tionality  factor  which  is  always  greater  than  zero;  s  is  called 
the  degree  of  reaction  or  reaction  coordinate.  It  indicates  the 
degree  to  which  the  reaction  has  taken  place  [19]. 

The  chemical  potential  (/x)  appearing  in  Eq.  (2)  (assuming 
a  pure  substance  in  a  reactive  mixture  of  ideal  gases)  [20]  is 
given  by 

/x  =  g°  RT  In  p  (4) 


At  equilibrium,  any  chemical  reaction  obeys  the  condition: 
dGr,p  =  0.  Consequently,  combining  Eqs.  (4)  and  (2),  the 
equilibrium  condition  can  be  expressed  in  terms  of  the  chem¬ 
ical  potential  as  follows: 


VCgc  +  VD#D  -  vA#a  -  yBgB  =  ~RT  ln 


( Pc)vc(pd)vd 
(Pa)Va(Pb)Vb 


The  term  in  the  left-hand  side  is  known  as  the  standard- state 
Gibbs  function  (A G°).  The  term  in  the  right-hand  side  is 
related  to  the  equilibrium  constant  ( Kp ) 

AG°  =  -RT\nKp  (6) 

The  standard- state  Gibbs  energy  corresponds  to  the  maxi¬ 
mum  work  that  can  be  drawn  from  a  fuel  cell  that  operates  at 
standard  conditions.  The  equilibrium  constant  can  be  found 
in  tables  or  usually  is  reported  in  terms  of  empirical  corre¬ 
lations.  The  partial  pressures  in  Eq.  (5)  can  be  expressed  in 
terms  of  their  equivalent  mole  fractions  (pi  =  Xi  P ),  thus  the 
equilibrium  constant  can  be  rewritten  as 

VQ  Vd 

Kp  =  lff%(p)vc+VD-VA-VB  (7) 

v  n  V  u 

AA  AB 

The  former  set  of  thermodynamic  relations  (Eqs.  (5)-(7))  are 
connected  with  the  fuel  cell  electrochemistry  by  the  reversible 
open  circuit  voltage  [21] 


Let  us  start  by  defining  the  following  chemical  reaction 

vaA  +  vbB  -<=>-  vcC  +  vdD  (1) 

where  the  capital  letters  indicate  the  chemical  constituents 
and  v  the  stoichiometric  coefficients  associated  with  the  bal¬ 
anced  reaction. 

Based  on  the  definition  of  Gibbs  function:  G  =  H  — 
TS ,  and  the  fundamental  thermodynamic  relation:  U  = 
U(S,  V,  Na,  Ab,  Nq ,  Ad),  the  Gibbs  energy  for  Eq.  (1)  is 


AG° 

~zF 


where  E°  is  the  maximum  voltage  obtainable  from  the  fuel 
cell,  z  the  number  of  electrons  transferred,  and  T  the  Faraday 
constant  {fF  =  9.6485  x  104Cmol_1). 

When  off-equilibrium  conditions  are  considered,  the  par¬ 
tial  pressures  at  equilibrium  conditions  (Eq.  (5))  must  be  re¬ 
placed  by  their  actual  values.  The  new  ratio  represents  the  off- 
equilibrium  concentrations  of  reactants  and  products  [22]. 
Thus  the  new  cell  potential  in  terms  of  actual  concentrations 
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(off-equilibrium)  is  given  by  the  well  known  Nernst  equation 
RT 

E  =  E° - In  K  (9) 

nT 

where  K  represents  the  ratio  of  products  to  reactants  at 
conditions  different  from  those  at  equilibrium. 

For  fuel  cells  operated  with  pure  hydrogen  (H2  + 
^02  =>  H2O),  the  Nernst  equation  becomes 


_  RTln  (  (pn2o)(P)  \ 
2T7  \(pu2)  Q Po2) ) 


(10) 


The  Nernst  equation,  however,  represents  an  idealized  situ¬ 
ation.  The  real  voltage  has  to  take  into  account  the  losses 
associated  with  activation  (ij ac),  ohmic  (77^),  and  concentra¬ 
tion  polarization  ( ijv )•  Moreover,  there  are  some  losses  oc¬ 
curring  in  the  fuel  cell  that  cannot  be  characterized  with  the 
above-mentioned  polarization  terms  (for  example,  leakage), 
but  in  general  their  contribution  is  insignificant.  Maloney  [23] 
studied  the  influence  of  different  rate-limiting  steps  in  fuel 
cell  performance.  Porous  gas  diffusion,  adsorption,  surface 
migration,  charge  transfer,  reaction  kinetics,  and  ohmic  resis¬ 
tance  were  considered.  Maloney  identified  ohmic  and  activa¬ 
tion  polarization  as  the  main  contributors  to  irreversibility. 
For  electrolyte- supported  cells,  the  main  polarization  is  at¬ 
tributed  to  the  ohmic  loss  [24].  For  anode-supported  cells, 
ohmic  polarization  is  small  because  of  the  thin  electrolyte 
used;  therefore,  the  main  contribution  comes  from  activa¬ 
tion  polarization.  Concentration  polarization  is  expected  to 
become  an  important  rate-limiting  step  at  high  current  den¬ 
sities  (>7500  Am-2  for  a  typical  membrane)  and  low-flow 
concentrations  (>80%  of  fuel  utilization)  [21]. 

Taking  into  account  the  three  main  sources  of  irreversibil¬ 
ity  occurring  in  the  fuel  cell,  the  cell  potential  is  given  by 


^Real  =  E  -  T)Ac,A  ~  0Ac,C  ~  ~  0V,A  ~  0V,C  (H) 

where  the  second  subscript  indicates  the  polarization  at  anode 
(A),  cathode(C),  and  electrolyte  (E). 

The  polarization  terms  are  discussed  next.  The  polariza¬ 
tion  model  assumes  that  the  reactions  occur  at  the  electrode¬ 
electrolyte  interface,  no  temperature  or  pressure  gradient 
throughout  the  electrode  thickness,  and  steady- state  condi¬ 
tions. 


2.7.  Activation  polarization 


where  j0  is  the  exchange  current  density,  a  the  electron  trans¬ 
fer  coefficient,  and  z  the  number  of  electrons  transferred  per 
reaction.  The  exchange  current  density  is  influenced  by  differ¬ 
ent  factors  which  are  not  clearly  understood  [12] .  Semiempir- 
ical  correlations  based  on  experimental  data  are  usually  used 
for  the  exchange  current  density;  however,  Campanari  et  al. 
[27]  showed  that  reported  values  in  the  literature  can  vary 
widely.  For  our  calculations  we  used  recommended  values 
by  Chan  et  al.,  but  we  recognize  that  more  work  needs  to  be 
done  on  this  area.  The  electron  transfer  coefficient  depends  on 
the  electrocatalytic  reaction  mechanism  and  typically  takes 
values  between  zero  and  one. 

The  Butler- Volmer  equation  is  solved  numerically  in 
this  paper.  The  most  commonly  used  method  for  solving 
non  linear  equations  is  the  Newton-Raphson.  The  MatLab® 
fzero  function  uses  the  zeroin  algorithm  for  solving  non¬ 
linear  equations.  The  zeroin  algorithm  combines  bisection, 
quadratic  interpolation,  and  secant  methods  for  speed  and 
reliability.  An  improved  version  of  this  algorithm  was  devel¬ 
oped  by  Moler  [28].  The  Mathematica®  FindRoot  function 
uses  Newton,  secant,  and  Brent’s  methods  for  finding  a  solu¬ 
tion  to  non-linear  equations.  Both  algorithms  were  used  for 
this  paper. 

At  high- activation  polarization,  the  second  term  in  the 
Butler- Volmer  equation  will  be  much  smaller  compared  to 
the  first  term  and  can  be  eliminated.  The  resulting  expression 
is  the  well  known  Tafel  expression 


0A 


Sin(t) 


(13) 


At  low-activation  polarization,  the  term  ( a  zTopf)/ RT  will 
be  much  less  than  unity  and  the  exponential  can  be  expanded 
as  a  Taylor  series 


’  " 

3' 

RT 

a 

(14) 


Neglecting  all  terms  with  an  order  higher  than  one,  and  solv¬ 
ing  for  rjA->  the  well  known  linear  current  potential  relation  is 
obtained 


0A 


(15) 


Activation  polarization  is  controlled  by  the  electrode  ki¬ 
netics  at  the  electrode  surface.  This  polarization  is  directly 
related  to  the  activation  barrier  that  must  be  overcome  by  the 
reacting  species  in  order  for  the  electrochemical  reaction  to 
occur.  The  electrode  reaction  rate  at  high  temperatures  (600- 
800  °C)  is  fast,  and  the  result  is  that  activation  polarization 
is  small,  which  represents  the  case  of  SOFCs.  Activation  po¬ 
larization  is  given  by  the  Butler- Volmer  equation  [25,26] 


exp 


et  rj Ac  zF 

RT 


—  exp 


Opl/Ac  zf\ 

RT  ) 


Whether  to  use  the  full  Butler- Volmer  equation  or  its  simpli¬ 
fied  forms  depends  on  the  accuracy  desired  and  the  expected 
activation  polarization.  Table  1  shows  the  ranges  at  which  the 
Tafel  and  the  linear  current  potential  are  within  a  5%  error 
relative  to  the  Butler- Volmer  equation.  These  ranges  were 
calculated  from  correlated  linear  expressions  (see  Eqs.  (16) 
and  (17)),  which  correspond  to  the  fitted  values  at  which 
the  ratio  j/j0  gives  an  error  of  5%  (Fig.  1)  compared  to 
the  Butler- Volmer  expression.  Table  1  is  in  good  agreement 
with  the  values  reported  by  Chan  et  al.  [11],  but  the  results 
presented  here  in  the  form  of  correlated  expression  gave  us 
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Table  1 


Activation  polarization  range  for  5%  accuracy  of  simplified  models  relative 
to  the  Butler- Volmer  equation 

P  (K) 

Tafel,  A 

Linear  current  potential,  r) a 

1073 

>0.2516 

<0.1006 

1173 

>0.2751 

<0.1100 

1273 

>0.2985 

<0.1194 

1373 

>0.3219 

<0.1287 

1473 

>0.3454 

<0.1382 

more  information  than  the  previous  results  presented  by  Chan 
et  al.  at  only  one  temperature. 

Tafel  error  correlation: 

r]A  >  0.0002345  T  -  8.06497  x  10“7  (16) 

Linear  current  potential: 

r]A  <  0.0000938  T  -  4.483069  x  10“7  (17) 


2.2.  Empirical  correlations  for  the  activation  loss 


It  is  common  to  use  semiempirical  correlations  to  calculate 
the  polarization  occurring  in  the  fuel  cell.  These  correlations 
are  preferred  because  of  their  simplicity;  however,  their  accu¬ 
racy  is  questionable  in  a  wide  range  of  operating  conditions. 
The  most  widely  used  correlations  are  the  ones  reported  by 
Achenbach  [15] 

Cathode: 


1 

~Rc 


4  T 

- rr 

RT 


(18) 


Anode: 


1 


^a,h2 


2T7  /  p\\, 

~RTrA’H2  l~/T 


1 

Ra.co 


2T7  /  pco 

—rh,CO  I  - 

RT  \  p 


m 


exp 


Ea 

RT 


(19) 

(20) 


where  pry,  =  0.21  p,  pry,  =  0.50p,  pro  =  0.01  p,  Er  = 
160kJ  mol-1,  Ea  =  1 10  kJ  mol-1,  rc  =  1.489  x  1010 


Temperature  (K) 


Table  2 


Resistance  correlations  for  calculating  the  activation  polarization;  R  = 
rexp(£y  T ) 


Reference 

Region 

r  (Q  m2) 

E(  K) 

Achenbach 

Cathode 

2.72  x  10“12 

19244 

Anode  (H2) 

3.07  x  lO"10 

13230 

Anode  (CO) 

5.82  x  10-10 

13230 

Hendriksen 

Cathode 

9.02  x  10“13 

20882 

Anode 

1.35  x  10“8 

8121 

Karoliussen 

Anode  and  cathode 

2.83  x  10“8 

8360 

Am  2,  L\,h2  —  2.128  x  108  Am  2,  and  rAco  = 

2.98  x  108 Am  2;  the  preexponential  factors  were 
calculated  under  the  assumption  of  R  =  0.1  £2  cm2.  The 
preexponential  factors  were  calculated  at  1000  °C;  therefore 
better  predictions  are  expected  near  this  temperature.  Achen¬ 
bach’ s  correlations  were  calculated  from  experimental  kinetic 
data  assuming  low  polarization.  A  quantitative  comparison 
between  three  different  empirical  correlations,  including 
Achenbach’s  expressions,  was  reported  by  Motloch  [22].  A 
summary  of  Motloch’ s  correlations  is  presented  in  Table  2. 
The  total  or  equivalent  resistance  is  calculated  considering 
the  fuel  cell  being  analogous  to  an  electrical  circuit,  where 
the  cathode  and  the  anode  are  connected  in  series.  If  carbon 
monoxide  and  hydrogen  are  considered  electrochemically 
active  fuels,  then  the  equivalent  resistance  is  calculated 
viewing  the  fuels  as  being  in  parallel:  R~l  =  Ra\i2  +  ^co- 
Motloch  found  that  the  variation  among  the  three  models 
is  minimal  at  high  temperatures,  T  >  1000  °C;  however,  at 
lower  temperatures,  T  <  800  °C,  these  expressions  produced 
unrealistic  results.  For  example  at  800  °C  (typical  value 
for  the  SOFC  operating  temperature),  the  overpotential  is 
r)  =  1 .0578  V,  assuming  a  current  density  of  7500  A  m-2. 

A  comparison  between  these  correlations  and  the  Butler- 
Volmer  equation  reveals  that  the  empirical  correlations  were 
reasonably  accurate  between  900  and  1473°C  (1173  and 
1273  K)  (Figs.  2  and  3).  At  higher  temperatures,  the  empirical 
expressions  gave  polarization  values  much  smaller  than  the 
Butler- Volmer  equation,  which  was  expected  based  on  the 
assumptions  used  for  these  correlations.  However,  at  lower 
and  higher  temperatures  (T  <  900  °C  and  T  >  1473  °C), 
the  numerical  values  reported  unrealistic  results  and  cannot 
be  used  for  any  practical  purpose.  Based  on  our  numerical 
findings,  we  recommend  the  use  of  the  full  Butler- Volmer 
equation  and  avoid  the  use  of  semiempirical  correlations, 
especially  at  low  temperatures.  Although  the  Tafel  and  the 
linear  current  potential  could  simplify  computational  time  in 
the  past,  today  the  implementation  of  algorithms  for  solving 
systems  of  algebraic  non-linear  equations  is  fast  enough, 
and  there  is  no  reason  to  sacrifice  accuracy  for  simplified 
relations. 

2.3.  Concentration  polarization 


Fig.  1.  Operating  region  for  Tafel  and  linear  current  potential  as  function  of 

temperature  and  activation  polarization  for  a  5%  error  relative  to  the  Butler-  Concentration  polarization  occurs  when  the  fuel  is  COn- 

Volmer  equation.  sumed  at  the  electrode-electrolyte  interface,  and  the  gas  con- 
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Fig.  2.  Comparison  between  the  total  activation  polarization  using  the 
Butler-Volmer  equation  (B)  and  empirical  correlations  reported  by  Achen- 
bach  [15]  (A),  and  Motloch  [22]  (H  =  Hendriksen,  K  =  Karoliussen);  cur¬ 
rent  density  j  =  3000  A  m-2  and  exchange  current  density  based  on  recom¬ 
mended  values  by  Chan  et  al.  [11]. 


centration  decreases  at  the  reaction  sites.  Concentration  po¬ 
larization  becomes  an  important  loss  at  high  current  densities 
and  small  fuel  concentrations  (>80%  of  fuel  utilization).  The 
main  factors  that  contribute  to  concentration  polarization  are 
diffusion  of  gases  through  the  porous  media  and  solution- 
dissolution  of  reactants  and  products.  The  transport  mecha¬ 
nisms  within  the  fuel  cell  are  governed  by  diffusion  transport 
and  Darcy’s  viscous  flow.  For  typical  operating  conditions  on 
an  SOFC,  diffusion  transport  will  dominate  and  convection 
transport  can  be  neglected  (Darcy’s  flow).  Molecular  diffu¬ 
sion  and  Knudsen  diffusion  describe  the  diffusion  transport 
through  the  porous  electrode.  Their  contribution  to  diffusion 
transport  is  tightly  related  to  the  microphysical  characteris¬ 
tics  of  the  porous  material  (i.e.,  porosity,  tortuosity,  pore  size, 
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Fig.  3.  Comparison  between  the  total  activation  polarization  using  the 
Butler-Volmer  equation  (B)  and  the  empirical  correlations  reported  by 
Achenbach  [15]  (A)  and  Motloch  [22]  (H  =  Hendriksen,  K  =  Karoliussen); 
current  density  j  =  7500  A  m-2  and  exchange  current  density  based  on  rec¬ 
ommended  values  by  Chan  et  al.  [11]. 


and  permeability).  The  diffusion  transport  in  a  porous  mate¬ 
rial  can  be  described  by  Fick’s  model,  the  dusty  gas  model  or 
the  Stefan-Maxwell  model.  Fick’s  model  is  used  more  fre¬ 
quently  because  it  is  simpler  to  implement  than  the  dusty  gas 
model  and  analytical  expressions  can  be  derived  more  easily 
[29] ;  the  Stefan-Maxwell  model  is  usually  discarded  because 
it  does  not  include  Knudsen  diffusion.  If  Knudsen  diffusion  is 
dominant,  the  dusty  gas  model  predictions  are  more  accurate 
than  those  from  Fick’s  model  [30].  The  reason  is  attributed 
to  the  way  the  overall  effective  diffusion  coefficients  are  de¬ 
termined  in  both  models. 


2.3.1.  Fick’s  model 

The  mass  transport  equation  of  a  single  component  fluid 
is  described  by 


g  Wji_P) 

RT  dt 


-  v  •  Nt  +  n 


(21) 


where  e  is  the  porosity,  A)  the  rate  of  mass  transport,  and  rz- 
the  reaction  rate  in  the  porous  electrode  [29] .  Assuming  the 
process  occurs  at  steady  state  and  that  the  electrochemical 
reaction  takes  place  at  the  anode-electrolyte  interface  the 
mass  transport  equation  becomes 


V  •  Ni  =  0 


(22) 


Fick’s  model  is  given  by 


1 

Ni  =  - 

RT 


3(y,-m 

3  z  ) 


(23) 


where  2\e  is  the  overall  effective  diffusion  coefficient  (see 
Appendix  A). 

For  the  binary  system,  H2-H2O,  the  rate  of  mass  transport 
as  given  by  Fick’s  model  is 


2\e  dCyHgjP) 
RT  dz 


(24) 


At  the  anode-electrolyte  interface  (z  =  l )  the  current  density 
produced  is  governed  by  the  rate  of  reactant  diffusing  into  the 
porous  anode  by 


(25) 


substituting  this  value  into  Eq.  (24),  the  following  boundary 
condition  is  obtained: 


RT 


2OTa,c 


J 


(26) 


By  solving  the  mass  transport  equation  for  the  system,  H2- 
H2O  (Eq.  (24))  with  boundary  conditions,  Eq.  (26)  and 
yn2  (z  =  0)  =  y{^lk  =  ypj  ,  the  composition  of  the  binary 
system  can  be  computed  at  the  electrode-electrolyte  inter¬ 
face 


2  /  Pa,e  P  dyn2  \  =  Q 
dz  v  RT  dz  ) 


(27) 
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Pu2  =/’h2-1ct  jz  (28) 

2  Da  e 

Noticing  that  pu2  +  pu2o  =  1 
T  RT 

Ph2o  =  Pu2o  +  2FDaJ  Z  (29) 

The  concentration  at  the  interface  is  therefore  calculated  in 
terms  of  the  partial  pressures  of  the  species  involved,  hydro¬ 
gen  in  this  case.  In  general,  for  any  system  containing  a  binary 
component  system  the  diffusion  transport,  according  to  Fick’s 
model,  is  given  by  Eq.  (28);  the  2  in  the  denominator  has  to  be 
modified  in  terms  of  the  number  of  electrons  transferred  when 
different  systems  are  considered.  The  Stefan-Maxwell  model 
is  applied  if  Knudsen  diffusion  is  assumed  to  be  insignificant. 
In  this  case,  the  same  procedure  developed  for  Fick’s  model 
is  used,  with  the  only  difference  being  that  the  overall  ef¬ 
fective  diffusion  coefficient  (2\e)  does  not  include  Knudsen 
diffusion.  The  analytical  treatment  of  multicomponent  sys¬ 
tems  using  Fick’s  model  becomes  excessively  complex  (from 
a  mathematical  point  of  view)  for  systems  greater  than  three 
components,  and  numerical  models  are  used  instead  [31]. 


2.3.2.  Dusty  gas  model 

The  dusty  gas  model  (DGM),  similar  to  Fick’s  model, 
includes  molecular  diffusion  and  Knudsen  diffusion  but  is 
based  on  the  Stefan-Maxwell  formulation.  The  rate  of  mass 
transport  according  to  the  DGM  is  given  by 


n 

+  E 

j=Ujfi 


yjNi  -  ytNj 


J^dy, 

RT  dz 


(30) 


The  DGM  assumes  that  the  pore  walls  consist  of  giant,  mo¬ 
tionless,  pseudo  molecules  (dust)  that  are  uniformly  dis¬ 
tributed  in  space.  The  flux  ratio  is  determined  using  the  Gra¬ 
ham’s  law  of  diffusion  in  gaseous  mixtures  whereas  in  Fick’s 
model  an  equimolar  counter  diffusion  concept  is  used  in¬ 
stead.  The  DGM  is  more  accurate  when  Knudsen  diffusion 
is  dominant,  because  at  these  conditions  the  assumption  of 
equimolar  counter  diffusion  is  no  longer  valid  [30]. 

Using  the  DGM,  the  compositions  for  the  system,  FU- 
H2O,  at  the  electrode-electrolyte  interface,  are  given  by 


yn2  +  Jh2o  =  1 


(31) 


d 2yu2  a 

dz2  ^h2— H2o,e 


l  -  a  yu2  ,  1 

- T - 

22h2— H20,e  £*H2 ,*,e  . 


with  the  following  boundary  conditions: 

yH2(0)  =  y\i2 


z= 0 


.  RT 

J2pF 


1  -ay\i2 

^H2-H20,e 


+ 


1 

-  ► 

At  2,&,e 


(32) 


(33) 

(34) 


where 


Mn2 

Mn2o 


The  derivation  for  ternary  systems  was  presented  by  Suwan- 
warangkul  et  al.  [30] .  In  this  paper,  the  DGM  for  multicompo¬ 
nent  systems  is  implemented  assuming  no  convection  trans¬ 
port  and  that  all  molar  fluxes  are  known  from  the  global  elec¬ 
trochemical  reactions.  Although  the  assumption  of  no  con¬ 
vection  transport  could  limit  the  usability  of  the  method  in 
some  cases,  it  considerably  reduces  its  complexity  for  com¬ 
putational  purposes.  For  instance,  the  method  proposed  by 
Zhu  and  Kee  [12]  requires  an  iterative  procedure  for  solving 
the  DGM,  while  the  derivation  of  analytical  expressions  be¬ 
comes  very  complicated  for  multicomponent  systems.  Using 
the  assumption  of  constant  pressure,  our  results  (multicompo¬ 
nent  systems)  are  quantitatively  comparable  to  those  reported 
by  others  [12],  while  the  problem  can  be  easily  treated  from 
a  computational  point  of  view. 


2.3.3.  Numerical  calculations  of  concentration 
polarization 

Concentration  polarization  is  defined  as  the  difference  be¬ 
tween  the  ideal  and  real  cell  voltage  (when  corrected  for 
ohmic  and  activation  polarizations).  The  ideal  voltage  is 
calculated  in  terms  of  the  bulk  concentration  in  the  stream 
channel,  whereas  the  real  cell  voltage  is  calculated  at  the 
electrode-electrolyte  interface.  For  the  binary  system,  H2- 
H2O,  the  concentration  polarization  is  given  by 


RT 

fir  = - In 

/c  1T 


/  yu2  •  yh20 

\  ya2o  ■  3'|  1, 


(35) 


where  the  superscript  I  indicates  inlet  conditions. 

Fick’s  model,  DGM  and  the  Stefan-Maxwell  model  were 
implemented  with  Mathemathica® , 1  and  the  results  are  pre¬ 
sented  next.  The  physical  parameters  of  the  porous  electrodes 
(i.e.,  porosity  and  tortuosity)  are  determined  from  recom¬ 
mended  values  by  other  authors  [17,32,33].  The  overall  dif¬ 
fusion  coefficient  is  calculated  first;  afterwards,  one  of  the 
three  models  is  implemented.  Finally,  the  resulting  flow  con¬ 
centration  (usually  in  terms  of  partial  pressures)  is  used  to¬ 
gether  with  the  definition  of  concentration  polarization. 

These  calculations  were  useful  to  identify  some  features 
of  the  system  behavior  to  diffusion  transport.  The  first  im¬ 
portant  feature  was  the  sensitivity  to  flow  concentration. 
The  shape  of  the  concentration  polarization  changed  from 
concave-down  to  concave-up  at  high-  and  low-flow  concen¬ 
tration,  respectively  (Figs.  4  and  5).  This  indicates  that  at 
low-flow  concentrations  the  diffusion  effects  became  more 
important  than  at  higher-flow  concentrations.  The  limiting 
current  at  50%  of  fuel  utilization  and  800  °C  was  found  at 


1  The  algorithms  used  for  these  calculations  are  available,  at  no  cost,  at: 

http://uweb.und.nodak.edu/~eduardo.hernandez.pacheco/Performance.htm. 
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Fig.  4.  Anode  concentration  polarization  for  an  anode-supported  cell  using 
Fick’s  (o),  dusty  gas  (•)  and  Stefan-Maxwell  (o)  models.  Anode  thickness 
<5  =  750  x  10-6  m,  T  =  800  °C. 

approximately  35  000  Am-2  (Fig.  5).  At  high-flow  concen¬ 
trations  the  DGM  and  the  Stefan-Maxwell  model  reported 
similar  results  (Fig.  4);  however  at  low-flow  concentrations, 
the  Stefan-Maxwell  model  gave  poor  predictions  with  re¬ 
spect  to  the  other  two  models  (Fig.  5).  If  methane  and  carbon 
monoxide  are  considered  electrochemically  active  fuels,  the 
diffusion  transport  will  play  an  important  role  in  determining 
the  performance  characteristics.  At  low-flow  concentration 
(80%  fuel  utilization),  the  diffusion  transport  of  methane  and 
carbon  monoxide  became  noticeably  lower  than  that  for  hy¬ 
drogen  (Fig.  6).  Lower  values  of  the  limiting  current  will 
be  obtained  with  heavier  compounds  because  diffusion  to 
the  reaction  sites  becomes  more  important.  This  indicates 
that  although  higher  power  densities  can  be  reached  with 
some  hydrocarbons,  concentration  polarization  will  be  an  im¬ 
portant  source  of  irreversibility.  These  numerical  results  re¬ 
vealed  diffusion  problems  with  fuel  cells  operated  with  heavy 
molecules. 
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Current  density  (A/sq.  m) 

Fig.  5.  Anode  concentration  polarization  for  an  anode-supported  cell  using 
Fick’s  (o),  dusty  gas  (•)  and  Stefan-Maxwell  (o)  models.  Anode  thickness 
<5  =  750  x  10"6  m,  T  =  800  °C. 
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Fig.  6.  Concentration  polarization  at  80%  of  fuel  utilization  considering  the 
direct  oxidation  of  hydrogen  ( — ),  methane  (o),  and  carbon  monoxide  (•). 
Anode  thickness  8  =  600  x  10-6  m ,  T  =  950  °C. 

Based  on  the  recommendations  given  by  Suwanwarangkul 
et  al.  and  the  results  shown  here,  we  concluded  that  the  DGM 
is  more  accurate  for  predicting  the  concentration  polarization 
and,  therefore,  was  the  method  chosen  for  our  performance 
calculations. 

2.4.  Ohmic  polarization 

Ohmic  resistance  polarization  occurs  in  the  electrode  ma¬ 
terials  (anode  and  cathode),  interconnects,  and  in  the  elec¬ 
trolyte.  This  loss  is  the  resistance  to  the  flow  of  electrons  in 
the  electrodes  and  ions  in  the  electrolyte.  The  ohmic  loss  is 
perhaps  the  major  loss  mechanism  in  an  electrolyte- supported 
fuel  cell  [24].  For  anode-supported  cells,  this  polarization 
will  not  have  a  big  effect  on  the  cell’s  performance.  The  elec¬ 
trode  ohmic  polarization  is  not  taken  into  account  because  of 
the  high  conductivity  of  the  electrodes;  therefore,  the  domi¬ 
nant  ohmic  loss  is  that  for  the  electrolyte.  Ohmic  polarization 
obeys  Ohm’s  law  and  is  given  by 

On  =  j  Rn  (36) 

where  Rq  represents  the  total  ionic  and  electronic  resis¬ 
tance  (expressed  in  terms  of  the  resistivity  of  each  mate¬ 
rial,  Rq  =  (p  L)  and  j  is  the  current  density.  The  resistivity 
of  the  different  materials  that  integrate  the  fuel  cells  is  re¬ 
quired  in  Eq.  (36).  A  summary  of  commonly  used  empirical 
correlations  for  determining  the  conductivity  and  resistivity 
of  Ni-YSZ  is  presented  in  Table  3  [22,34,35]. 


Table  3 

Resistivity  ip  (£2  cm)  and  electrical  conductivity  a  (£2_1  cm-1)  parameters 
for  YSZ 


A 

B 

Equation 

7.86  x  10 S/T 

-10556 

A  +  exp  (B/T) 

190 

-9281 

A  +  exp  (B/T) 

334 

-10300 

A  +  exp  (B/T) 

0.3685 

10300 

(A  +0.002838  exp 

0.00294 

10350 

(AexpCB/T))-1 
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These  correlations  showed  an  almost  perfect  agreement  at 
medium  temperatures  (^600  °C)  and  a  small  discrepancy  at 
high  temperatures  (>900  °C).  The  variation  of  these  correla¬ 
tions  could  be  attributed  to  the  different  composition  used  in 
the  electrolyte.  We  have  chosen  the  correlation  reported  by 
Bessette  because  of  its  consistency  with  the  other  correlations 
in  the  entire  range  analyzed  and  its  simplicity. 


3.  Results 

The  performance  of  the  fuel  cell  was  determined  using  the 
macro-level  model.  For  pure  electronic  conductors  as  elec¬ 
trodes,  it  is  logical  to  assume  that  the  reactions  take  place 
at  the  electrode-electrolyte  interface.  However,  for  com¬ 
posite  electrodes  as  Ni-YSZ  (cermet),  a  micro-level  model 
will  be  required.  Although  the  assumption  of  a  pure  elec¬ 
tronic  conductor  does  not  apply  to  cermet  electrodes,  the 
numerical  results  reported  by  Chan  et  al.  [11]  showed  that 
for  anode- supported  cells  with  thick  anodes  (>600  pan),  the 
main  contribution  to  polarization  occurred  near  the  interface. 
The  SOFC  membranes  modeled  in  this  paper  have  anodes 
of  600-750  |mm.  This  anode  thickness  corresponds  to  values 
found  in  commercial  membranes  that  are  being  investigated 
at  the  EERC  and  UND  for  the  BG-SOFC  system.  Porosity 
and  tortuosity  can  be  determined  experimentally;  however, 
recommended  values  were  used  for  this  model.  Corrected 
values  of  these  parameters  will  be  determined  from  perfor¬ 
mance  experiments  in  future  work;  these  parameters  can  be 
corrected  in  order  to  fit  the  performance  curves  obtained  ex¬ 
perimentally.  The  algorithm  to  compute  the  polarization  in 
the  SOFC  proceeds  as  follows.  The  activation  polarization 
was  calculated  using  the  Butler- Volmer  equation.  The  zeroin 
algorithm  was  used  in  this  paper,  the  exchange  current  den¬ 
sity  was  based  on  values  recommended  by  Chat  et  al.  and 
Zhu  and  Kee.  The  DGM  was  implemented  for  determining 
the  diffusion  transport  in  the  anode.  The  numerical  values  for 
the  cathode  electrode  were  found  insignificant  compared  to 
those  at  the  anode.  A  sample  code  is  presented  in  Appendix  B 
to  determine  the  concentration  polarization  at  the  interface  of 
a  binary  system.  The  dusty  gas  model  requires  the  knowledge 
of  the  fuel  composition  at  the  bulk  stream.  This  composition 
was  determined  based  on  the  fuel  utilization  desired.  The 
ohmic  polarization  resistance  was  calculated  using  a  semiem- 
pirical  correlation  reported  by  Bessette.  At  the  electrodes,  an 
insignificant  resistance  is  expected  because  of  the  high  con¬ 
ductor  material.  The  algorithm  permitted  us  to  investigate  the 
optimal  operating  conditions  for  optimal  performance. 

The  macro-level  model  was  implemented  on  a  membrane 
operated  at  80%  of  fuel  utilization  and  800  °C.  The  fuel  is 
assumed  to  be  pure  hydrogen.  The  maximum  power  was 
obtained  at  approximately  3600  Wm-2  and  a  current  den¬ 
sity  of  9600  Am-2  (Fig.  7).  At  this  value,  the  polarization 
contribution  was  distributed  as  follows:  (1)  cathode  activa¬ 
tion  polarization  (0.2976  V),  (2)  anode  activation  polarization 
(0.1503  V),  (3)  ohmic  polarization  (0.0436  V)  and  (4)  anode 
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Fig.  7.  Performance  characteristics  for  hydrogen  at  80%  of  fuel  utilization 
and  800  °C.  Cathode  activation  polarization  (CA),  anode  activation  polar¬ 
ization  (AC),  anode  concentration  polarization  (AD)  and  electrolyte-ohmic 
polarization  (OP). 


concentration  polarization  (0.0385  V).  The  cathode  activa¬ 
tion  polarization  reported  the  highest  values  because  of  the 
low  exchange  current  density  used  (j0  =  200 Am-2).  Al¬ 
though  anode  concentration  polarization  did  not  represent  an 
important  loss  at  this  current  density,  at  higher  values,  this 
polarization  increased  very  rapidly,  especially  near  the  limit¬ 
ing  current  density  of  1 8  000  A  m-2 .  Fuel  utilization  is  tightly 
related  to  the  concentration  polarization,  because  more  fuel 
is  demanded  at  high  fuel  utilization,  more  concentration  po¬ 
larization  is  expected.  Fuel  cells  usually  operate  at  70-80% 
of  fuel  utilization  for  a  compromise  between  performance 
and  operational  costs.  Because  of  the  thin  electrolyte  used 
(10  pan),  ohmic  polarization  had  the  smallest  contribution  to 
the  voltage  loss,  but  for  different  types  of  cells  (electrolyte- 
supported  cells),  this  polarization  represents  the  main  factor 
for  reducing  the  performance  of  the  cell. 

Temperature  also  plays  an  important  role  in  determining 
the  performance  characteristics  of  the  fuel  cell.  As  Fig.  8 
shows,  the  maximum  power  density  shifted  to  a  lower  value 
when  the  temperature  decreased.  The  performance-jump 
from  800  to  700  °C  was  smaller  than  the  jump  from  700 
to  600  °C.  The  main  reason  for  this  is  because  at  low  tem¬ 
peratures,  both  ohmic  polarization  and  concentration  polar¬ 
ization  reported  higher  values  of  voltage  loss  (Fig.  9(a)  and 
(b)).  This  is  especially  notorious  for  ohmic  polarization.  At 
low  values  of  temperature  and  high  values  of  current  density, 
the  ohmic  polarization  increased  more  rapidly  than  the  other 
losses.  Consequently,  the  change  in  power  density  at  600  °C 
had  a  bigger  impact  than  the  change  at  higher  temperatures 
(700-800  °C).  At  high  temperatures  (>700  °C)  the  main  con¬ 
tribution  to  voltage  loss  came  from  activation  polarization 
(Fig.  9(c)  and  (d)).  As  temperature  and  current  density  in¬ 
creased,  the  activation  polarization  also  increased  but  not  as 
drastically  as  the  ohmic  polarization  did  at  low  temperatures 
(^600  °C). 
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Fig.  8.  Performance  characteristics  for  hydrogen  at  80%  of  fuel  utiliza¬ 
tion  and  different  temperatures.  (A)  T  =  600°C,  (B)  T  =  700°C,  (C) 
T  =  800  °C. 


Fig.  10.  Performance  characteristics  of  three  electrochemically  active  fuels: 
hydrogen,  methane,  and  carbon  monoxide.  The  numerical  calculations  were 
made  at  80%  of  fuel  utilization. 


The  effect  of  using  different  electrochemically  active  fuels 
is  investigated  next.  Although  the  direct  oxidation  of  hydro¬ 
carbons  is  not  possible  to  date  with  conventional  Ni-YSZ  an¬ 
odes,  there  are  many  research  groups  developing  new  anode 
materials  because  of  the  abundant  benefits  of  this  technology 
[36].  As  mentioned  before,  the  gases  leaving  the  gasifier  will 
be  composed  of  hydrogen,  methane,  carbon  monoxide,  car¬ 
bon  dioxide,  and  small  amount  of  tars  (1-5%)  [3].  Methane 
is  more  likely  to  reform  into  hydrogen  and  carbon  monoxide 
rather  than  being  directly  oxidized.  Carbon  monoxide  can  be 
directly  oxidized  on  Ni-YSZ  anodes  but  the  water-gas  shift 
reaction  is  more  favorable  in  the  presence  of  water.  The  de¬ 
tailed  understanding  of  all  the  different  reaction  mechanisms 
that  occur  in  the  fuel  cell  is  fundamental  for  a  complete  analy¬ 


sis  of  the  performance  characteristics  of  a  fuel  cell.  However, 
for  a  generic  analysis,  some  simplifications  can  be  made.  For 
instance,  for  a  H2-CO  system,  the  water-gas  shift  reaction 
is  assumed  to  be  always  at  equilibrium  and  hydrogen  can  be 
considered  the  only  fuel  in  the  system  [37].  For  a  more  com¬ 
plex  system  like  natural  gas  (85%  methane),  the  methane 
can  be  assumed  to  be  internally  reformed  at  high  temper¬ 
atures  (500-800  °C),  and  the  reaction  rates  will  determine 
the  amount  of  hydrogen  being  produced  in  the  system  [38]. 
Fig.  10  shows  the  performance  obtained  from  a  hypothetical 
fuel  cell  operated  with  pure  methane  and  pure  carbon  monox¬ 
ide.  The  benefits  obtained  from  a  higher  Nernst  potential  in 
the  case  of  methane  are  rapidly  overshadowed  by  diffusion 


Fig.  9.  Polarization  curves  ( V  [V])  as  a  function  of  temperature  (T  [K])  and  current  density  (j  [Am~2])  for  hydrogen  at  80%  of  fuel  utilization,  (a)  Ohmic 
polarization,  (b)  concentration  polarization,  (c)  anode  activation  polarization,  and  (d)  cathode  activation  polarization.  Current  density  =  0-18 000  Am-2, 
temperature  =  600-800  °C  (800-1073  K). 
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problems  (concentration  polarization).  The  window  of  op¬ 
eration  for  this  hypothetical  cell  will  be  small.  Although  the 
benefits  of  this  technology  (direct  oxidation  of  hydrocarbons) 
are  vast,  new  problems  will  arise.  For  now,  internal  reforming 
and  partial  oxidation  will  be  a  more  feasible  technology  for 
the  use  of  hydrocarbons  in  SOFCs  [39]. 


4.  Final  remarks 

A  macro-level  model  was  used  to  determine  the  per¬ 
formance  characteristics  of  an  anode-supported  cell.  The 
Butler- Volmer  equation  was  recommended  for  the  activation 
polarization;  however,  for  simplified  models  (Tafel  and  linear 
current  potential),  ranges  of  5%  error  relative  to  the  Butler- 
Volmer  equation  were  reported  in  the  form  of  two  correlated 
expressions.  Concentration  polarization  was  calculated  using 
the  DGM.  An  algorithm  of  this  model  implemented  on  a 
binary  component  system  is  included  in  this  paper.  For  mul¬ 
ticomponent  systems  the  algorithms  are  available  from  the 
authors.  Cathode  concentration  polarization  was  not  included 
because  the  diffusion  transport  contribution  was  insignificant 
because  of  the  thin  cathode  used  (40  [xm).  For  the  ohmic  po¬ 
larization,  different  empirical  correlations  showed  minimal 
variation,  especially  at  SOFC  operating  temperatures  (600- 
800  °C).  The  contribution  to  voltage  loss  was  analyzed  for 
hydrogen  at  80%  of  fuel  utilization  and  800  °C  on  a  anode- 
supported  cell.  A  surface  chart  for  each  polarization  term  as 
functions  of  temperature  and  current  density  was  presented. 
The  model  was  also  implemented  for  carbon  monoxide  and 
methane,  assuming  they  were  electrochemically  active  fuels. 
The  results  showed  that  concentration  polarization  will 
be  an  important  problem  associated  with  this  technology. 
Our  objective  was  to  provide  a  set  of  simple  computational 
tools  for  calculating  the  performance  characteristics  of  an 
SOFC  and  determining  the  operating  conditions  for  optimal 
performance.  Fuel  cells  usually  operate  to  the  left  of  the 
maximum  power  density,  because  of  the  better  voltage 
efficiency  and  stability  of  the  membrane.  However,  at  low 
current  densities  (left  of  the  peak  power  density)  carbon 
deposition  is  expected  (for  the  case  of  dry  methane  or  low 
water  content),  which  will  limit  the  window  of  operating 
conditions.  This  analysis  and  the  tolerance  of  SOFCs  to  tars 
is  being  investigated  at  the  EERC  and  will  be  presented  in  the 
future.  Our  results  showed  the  operating  conditions  at  which 
the  optimal  performance  is  obtained.  However,  in  reality, 
the  optimal  performance  is  a  compromise  among  maximum 
power  density,  operational  stability,  and  capital  costs. 
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Appendix  A.  Diffusion  coefficients 


Diffusion  in  porous  media  is  usually  described  by  molecu¬ 
lar  diffusion  or  Knudsen  diffusion.  Knudsen  diffusion  occurs 
when  the  diameter  pores  are  small  compared  with  the  mean 
free  path  of  the  gas  molecules.  Ordinary  diffusion  occurs 
when  the  pore  diameter  is  large  compared  with  the  mean  free 
path  of  the  gas  molecules.  The  mean  free  path  represents  an 
average  of  the  distance  that  the  molecules  can  travel  before  a 
collision  occurs.  The  Knudsen  number  Kn  =  X/d  (where  X 
is  the  mean  free  path  length  and  d  the  pore  diameter)  charac¬ 
terizes  the  diffusion  process.  For  Kn  <<C  1  ordinary  diffusion 
dominates,  for  Kn  1  Knudsen  diffusion  dominates.  For 
SOFCs,  both  Knudsen  and  ordinary  diffusion  processes  have 
to  be  considered  since,  in  general,  Kn  ~  1 . 

For  straight  and  round  pores  [11],  Knudsen  diffusion  is 
given  by 


2  /  8000 R  T 

3ry  Mi 


(A.l) 


2s 

r  =  - 

SaPb 


(A.  2) 


where  Sa  is  the  surface  area  of  the  porous  solid  (m2  kg  - 1 ),  pb 
the  bulk  density  of  the  solid  particle  (kg  m-3),  s  the  porosity 
material,  and  M  the  molecular  mass  (kg  kmol-1). 

In  order  to  account  for  the  tortuosity  of  the  material,  the 
Knudsen  coefficient  has  to  be  modified  in  terms  of  an  effective 
coefficient  [11,12]: 


(A. 3) 


4  £  /8 RT 
3r§y  tv  Mi 


(A.4) 


where  e/§  represents  the  porosity :tortuosity  ratio.  Zhang  et 
al.  [32]  proposed  a  new  way  to  predict  tortuosity  in  catalyst 
pellets  based  on  fractal  geometries;  a  hard  spheres  model  is 
often  used  for  determining  tortuosity.  The  accurate  determi¬ 
nation  of  this  parameter  is  out  of  the  scope  of  this  paper  and 
a  recommended  values  of  2. 0-6.0  were  used  instead. 

The  binary  diffusion  coefficients  can  be  determined  from 
the  Chapman-Enskog  theory  [40] 


3  J4: VkfJWj 

16  rut  erf 
6  / 


(A. 5) 


If  gases  are  assumed  to  be  ideal,  the  former  equation  can  be 
simplified  to 


0.0026  T3/2 
P  Mlf  fj 


(A.  6) 


where  /d  is  assumed  to  take  a  value  of  one,  and  n  is  deter¬ 
mined  using  the  ideal  gas  law.  When  using  these  equations, 
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Table  A.  1 


Lennard-Jones  potentials 


n2 

O2 

ch4 

h2o 

CO 

h2 

co2 

<*i 

3.798 

3.467 

3.758 

2.641 

3.690 

2.827 

3.941 

£i/k 

71.400 

106.700 

148.600 

809.100 

91.700 

59.700 

195.200 

Table  A.2 


Collision  integral  constants 


A 

B 

C 

D 

E 

F 

G 

H 

1.06036 

0.15610 

0.19300 

0.47635 

1.03587 

1.52996 

1.76474 

3.89411 

the  following  observations  have  to  be  considered: 


Mj  j  =  2  x 


1  1  'i 

Mi  +  Mj  ) 


°i  +  ai 
2 


(A. 7) 
(A. 8) 


AC  EG 

rB  +  exp(Z)r)  +  exp(Fr)  +  Ht 


X  = 


kT 


1/2 


(A. 9) 

(A.  10) 
(A.11) 


where  k  =  1.38066  x  10  23  (JK  l)  is  the  Boltzmann’s  con- 

o 

stant,  o' /  j  (A)  the  characteristic  length,  Sfj  (K)  the  char¬ 
acteristic  Lennard-Jones  length,  and  the  collision  inte¬ 
gral  based  on  the  Lennard-Jones  12-6  potential.  Values  for 
the  characteristic  lengths  are  reported  in  Table  A.l  and  the 
constants  appearing  in  the  collision  integral  are  reported  in 
Table  A.2  (taken  from  Reid  et  al.  [40]). 

Similar  to  the  effective  Knudsen  coefficient,  the  binary 
diffusion  coefficient  has  to  be  modified  in  order  to  account 
for  the  tortuosity  of  the  material, 


(A.  12) 


Because  ordinary  and  Knudsen  diffusion  might  occur  simul¬ 
taneously  an  overall  effective  diffusion  coefficient  is  given 


(A.  13) 


Todd  et  al.  [41]  compared  different  approximations  for  the 
determination  of  molecular  diffusion  for  binary  systems.  The 
approximations  were  compared  with  the  scarce  experimental 
data  available.  He  recommended  the  Fuller  et  al.  [40]  method 
because  it  was  consistently  accurate  over  the  entire  temper¬ 
ature  range  analyzed.  Based  on  this  method,  Table  A. 3  con¬ 
tains  a  set  of  rational  approximations  for  the  binary  diffusion 
on  a  system  that  contains  H2,  H2O,  CO,  CO2  and  CH4;  the 
method  of  rational  approximation  has  been  used  previously 
for  the  determination  of  thermophysical  properties  of  gases 
[42,43]. 


ao  +  a\T  +  C12T2 
1+biT 


Appendix  B.  Dusty  gas  model  (sample  code) 

The  dusty  gas  model  (DGM)  for  a  binary  component  sys¬ 
tem  (H2-H2O)  was  implemented  using  Mathematica®.  The 
fuel  utilization  is  set  to  80%  for  three  current  densities  (3000, 
7500,  and  15  000  Am-2)  and  a  thickness  anode  of  600  pm. 
The  diffusion  coefficients  are  calculated  using  the  method 
proposed  by  Fuller  et  al.  [40].  The  concentration  polarization 
is  calculated  once  the  diffusion  values  are  known  (last  line  of 
code). 


Table  A.3 


Binary  diffusion  coefficients  as  function  of  temperature 


a0 

a\ 

<32 

b\ 

R2 

SSE 

h2-h2o 

-0.152275 

0.001572 

7.031465e-06 

0.000109 

0.9999 

0.001755 

h2-co 

-0.131687 

0.001363 

6.114805e-06 

0.000109 

0.9999 

0.001515 

h2-co2 

-0.106085 

0.001127 

5.183760e-06 

0.000112 

0.9999 

0.001226 

h2-ch4 

-0.113444 

0.001201 

5.509623e-06 

0.000111 

0.9999 

0.001306 

h2o-co 

-0.041895 

0.000443 

2.034715e-06 

0.000111 

0.9999 

0.000482 

h2o-co2 

-0.034615 

0.000361 

1.627449e-06 

0.000110 

0.9999 

0.000396 

h2o-ch4 

-0.042982 

0.000450 

2.045609e-06 

0.000110 

0.9999 

0.000492 

co-co2 

-0.026500 

0.000280 

1.281830e-06 

0.000111 

0.9999 

0.000304 

CO-CTL 

-0.035039 

0.000370 

1.696197e-06 

0.000111 

0.9999 

0.000403 

co2-ch4 

-0.029049 

0.000305 

1.387365e-06 

0.000111 

0.9999 

0.000332 
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Constants  and  operating  conditions 
R  =  8.314;  F  =  96485;  i  =  0.00075;  cf  =  101325;  p  =  1; 
yAin  =  0.2;  yBin  =  0.8;  T  =  1073; 

Binary  diffusion:  H2-H2O  system 
Fuller  method: 

MA  =  2.016;  VA  =  6.12; 

MB  =  18.015;  VB  =13.10; 

VAB  =  (VA1/3 4  +  VB1/3)2; 

MAB  =  2(MA-1  +  MB-1)-1; 
a  =  1— Sqrt[MA/MB]; 

£  =  0.30;  §  =  6.0; 

DAB  =  (e/f)((0. 000143  TL75)/(p  MAB  5 6 7 8 9  VAB)); 

DAB  =  DAB/ 100 
Knudsen  diffusion: 
r  =  0.0000005; 

DAK  =  (e/§)  (97.0  r  Sqrt[T/MA]) 
j  =  {3000,7500,15000};  In  =  Length{j}; 
dusty  gas  model: 

Flatten [Table[  {y  1  =Evaluate[  y[l]  /. 

NDSolvef  |y”[x]  +  (a/DAB)  ((1  —  a  y[x])/ 

DAB  +  1/DAK)-1  y’[x]2  ==  0, 
y  [0]  =  yAin, 

y’[0]  =  —  <j[[m]]  R  T)/(2  p*cf  F)  ((1  -  a  yAin)/(DAB) 
+  1/DAK)-1}, 
y,  {x,0,£}]];  y2  =1  -  yl; 

-  ((R  T)/(2  F))  Log[(yl  yBin)/(y2  yAin)]},{m,l,ln}]]; 
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